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The  purpose  of  this  study  was  to  examine  the  feasibility 


of  modeling  damage  initiation  in  composite  materials  using 


finite  elements.  Composites  have  become  an  important  part 


of  aircraft  design.  The  ability  to  model  and  predict  damage 


initiation  will  help  in  the  design  as  well  as  the  failure 


tolerance  analysis  of  structures  comprised  of  composite  ma¬ 


terials.  An  existing  finite  element  code,  NOSAPM,  was  modi¬ 


fied  to  include  a  failure  criteria  and  stiffness  reduction 


techniques.  Along  with  the  headaches  of  changing  a  compli¬ 


cated  code  such  as  NOSAPM  were  the  problems  of  converting 


the  mounds  of  generated  data  into  something  useful.  The 


post-processing  routines  originally  written  for  NOSAPM  are 


either  extinct  or  well  hidden  within  the  bowels  of  the 


Wr ight-Patterson  AFB  computer  network. 
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Abstract 

The  purpose  of  this  study  was  to  analyze  damage  initia¬ 
tion  in  the  short  beam  shear  test  for  composite  materials. 
This  thesis  proposes  a  means  of  predicting  damage  initiation 
in  a  composite  beam  using  finite  element  techniques  and  a 
continuum  mechanics  failure  criteria.  The  first  objective 
was  to  accurately  determine  the  stress  field.  The  second 
objective  was  to  employ  a  damage  criteria  based  on  the  ele¬ 
mental  stress  field.  The  final  objective  was  to  evaluate 
the  effectiveness  of  using  an  elemental  stiffness  reduction 
to  study  damage  initiation  and  progression.  The  finite  ele¬ 
ment  code  used  was  NOSAPM. 

The  complex  stress  fields  associated  with  this  problem 
were  accurately  obtained  with  the  exception  of  the  shear 
stresses.  The  high  shear  gradients  required  the  use  of  a 
very  fine  mesh  which  was  not  possible  with  the  computer 
resources  available.  The  Tsai-Hill  failure  criteria  was 
successfully  Incorporated  into  NOSAPM.  Also,  an  elemental 
stiffness  routine  was  derived  and  verified  with  a  single 
element  compression  case.  With  the  beam  test  model,  the 
stiffness  reduction  technique  did  not  dramatically  affect 


the  material  with  fibers  oriented  ninety  degrees.  However, 
the  zero  degree  fiber  orientation  case  resulted  in  stress 
contour  changes  which  suggested  a  spreading  of  the  load 
through  an  element.  This  spreading  only  occurred  with  a 
stiffness  reduction  greater  than  ninety-five  percent.  At 
high  loads  and  stiffness  reductions,  the  model  breaks  down 
as  the  numerical  solutions  result  in  element  instabilities. 

The  most  significant  recommendations  are  to  refine  the 
mesh  using  a  two  dimensional  code  and  to  reduce  the  stiff¬ 
ness  only  in  the  area  best  representing  the  failure  mode. 
Also,  load  shifting  should  be  investigated  to  prevent  the 
elemental  instabilities.  Finally,  a  study  using  an  lndentor 
model  instead  of  a  point  load  should  be  accomplished. 


Damage  Initiation  in  the 
Short  Beam  Shear  Test  of  Composite  Materials 
using  Stiffness  Reduction 


Li.  Introduction 

An  important  topic  in  the  study  of  the  mechanics  of 
materials  is  the  characterization  of  damage.  The  ability  to 
predict  damage  enables  an  engineer  to  account  for  it  and  to 
make  the  required  design  changes  as  appropriate.  Conven¬ 
tional  materials  such  as  steels  and  metal  alloys  have  fairly 
well  established  damage  behaviors.  Composite  materials, 
however,  often  do  not  fail  in  the  same  mode  as  their  metal 
counterparts.  With  the  increased  use  of  composites  in  air¬ 
craft  design,  better  techniques  are  needed  to  evaluate  and 
predict  damage  to  such  composite  materials  as  graphite/epoxy 
and  the  more  exotic  metal  matrix  materials.  This  thesis  is 
an  attempt  to  understand  one  aspect  of  damage  modeling,  the 
initiation  of  damage  caused  by  a  bending  load. 


Material  properties  are  determined  by  testing.  Pure 
tension  and  compression  tests  have  been  used  predominantly 


due  to  the  simplicity  of  specimen  loading  and  geometry. 

These  tests  provide  an  understanding  of  failure  by  reducing 
the  number  of  variables  considered.  However,  biaxial  tests 
like  the  short  beam  shear  test  provide  a  more  complex  stress 
field  and  require  more  insight  into  the  modes  of  material 
damage.  This  short  beam  test  is  now  being  scrutinized  to 
determine  its  effectiveness  as  a  means  of  measuring  inter¬ 
laminar  shear  strength.  The  complex  stress  field  of  the 
test  has  produced  local  failures  that  occur  prior  to  the 
expected  shear  delaminations  thus  casting  suspicion  on  test 
results.  Comprehensive  modeling  of  damage  initiation  and 
progression  is  essential  in  developing  an  understanding  of 
this  phenomena. 

In  order  to  make  an  assessment  of  composite  damage  and 
how  to  model  it  one  must  first  understand  how  it  originates. 
To  date,  test  observations  are  inconclusive  as  catastrophic 
failure  occurs  too  rapidly  to  allow  identification  and  mea¬ 
surement  of  local  and  global  damage.  Also,  the  microscopic 
details  of  the  physics  of  composite  damage  have  not  been 
fully  explained  or  modeled.  This  thesis  proposes  a  means  of 
predicting  damage  initiation  in  a  composite  beam  using  fi¬ 
nite  element  techniques  and  a  continuum  mechanics  failure 
criteria.  The  first  objective  is  to  accurately  determine 
the  stress  field  for  the  short  beam  shear  test  of  a  compos- 


ite  beam.  A  convergence  study  of  the  required  mesh  will  be 
conducted  using  the  finite  element  code  NOSAPM.  The  second 
objective  is  to  employ  a  continuum  mechanics  damage  criteria 
that  predicts  damage  based  upon  the  elemental  stress  field. 
The  third  objective  is  to  determine  the  effectiveness  of 
using  stiffness  reduction  to  study  damage  initiation  and 
progression.  Stiffness  reduction  methods  have  been  employed 
successfully  in  simple  stress  fields.  It  is  of  fundamental 
interest  to  determine  if  these  methods  can  also  be  success¬ 
fully  applied  to  complex  stress  fields  to  accurately  charac¬ 
terize  the  experimentally  observed  local  and  global  failure 
modes.  The  results  of  this  thesis  will  serve  as  a  basis  to 
determine  the  most  fruitful  approach  that  should  be  pursued 
by  follow-on  efforts  to  characterize  damage  in  composite 
mater ials  . 
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The  question  of  how  to  model  damage  in  composite  materi¬ 
als  is  not  a  new  one.  Various  methods  have  been  employed  to 
try  to  describe  the  phenomena  of  material  failure.  To  begin 
with,  failure  must  first  be  defined.  The  point  when  a  load 
exceeds  the  capacity  of  a  structure  to  elastically  support 
it,  normally  called  yielding,  is  the  most  widely  applied 
failure  criteria.  However,  for  a  complex  loading  a  more 
convenient  measure  would  be  a  relation  of  some  failure  cri¬ 
teria  to  the  entire  stress  field  of  the  structure.  This 
stress  field  could  be  most  easily  defined  by  using  simple 
beam  theory.  But,  beam  theory  breaks  down  with  large  dis¬ 
placements  and  also  fails  to  account  for  localized  loading 
effects.  The  strain-displacement  and  stress-strain  rela¬ 
tionships  of  finite  element  techniques,  however,  do  provide 
a  way  to  overcome  these  beam  theory  limitations.  Of  course, 
finite  elements  do  depend  significantly  on  the  modeling  of 
the  structure  and  the  types  of  elements  used. 

Alternatively,  one  can  solve  the  elasticity  contact 
problem  for  the  complex  stress  fields  produced  by  a  rigid 
indenter.  Such  a  mixed  boundary  value  problem  requires  the 
application  of  integral  transforms.  While  these  techniques 


mm 


may  be  more  exact,  integral  transform  solutions  for  multi¬ 
layered  structures  can  not  be  obtained  in  closed  form  thus 
prompting  one  to  use  finite  elements.  Regardless,  all  nu¬ 
merical  techniques  must  be  validated  by  physical  testing,  if 
possible,  to  determine  accuracies  and  limitations  of  the 
model ing . 

Damage  must  be  defined  using  some  criteria.  Exceeding 
this  criteria  indicates  damage  will  occur.  As  with  conven¬ 
tional  materials,  there  exist  many  theories  to  define  damage 
in  composite  materials.  The  criteria  used  for  this  research 
is  based  upon  the  longitudinal,  transverse,  and  shear 
strengths  of  the  material.  It  was  chosen  because  it  uti¬ 
lizes  the  three  dimensional  state  of  stress  throughout  the 
body.  For  conventional  materials  this  strength  is  best  de¬ 
scribed  by  yield  criteria  such  as  the  Tresca  (Maximum  Shear 
Stress)  and  Von  Mises  ( Distort ional  Energy)  criteria 
(10:134,139).  By  using  the  element's  stress  fields  with  a 


combination  of  material  strengths  one  obtains  an  indication 
of  potential  failure. 

Methods  to  determine  the  strength  of  a  composite  materi¬ 
al  are  shown  by  both  Jones  (8:71,83)  and  Tsai  (17:277,325). 
Modern  theories  include  the  maximum  stress  theory,  the  maxi¬ 
mum  strain  theory,  the  Tsai-Wu  tensor  theory,  and  the 
Tsal-Hlll  theory.  Each  of  these  theories  have  strengths  and 
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weaknesses,  so  application  of  any  must  be  carefully  consid¬ 
ered.  A  brief  look  at  each  of  these  theories  is  necessary 
to  establish  the  one  best  suited  for  this  research. 

According  to  Jones,  the  maximum  stress  theory  requires 
that  the  stresses  in  principal  material  directions  be  less 
than  material  strengths  in  those  directions  (8:72).  Thus, 

Xc  <  <7j  <  X, 

yc<°2<Y, 

Ze<a3<Zt 

|tI2|<s  (l) 


b-: 


i 


j 

'j 


»,i 

v.v 


where  X,  Y,  and  Z  indicate  material  strengths  in  the  respec¬ 
tive  directions,  S  is  the  shear  strength,  sigma  is  the 
stress  in  the  material  principal  directions,  and  the  sub¬ 
scripts  "t"  and  "c"  represent  tension  and  compression,  re¬ 
spectively.  In  a  similar  manner  the  maximum  strain  theory 
compares  the  strains  in  the  principal  material  directions 
with  the  maximum  allowable  strains  in  those  directions. 
However,  these  theories  do  not  take  into  account  any  com¬ 
bined  effects  of  the  complex  stress  fields. 

Wu  (21:361-400)  published  a  comparison  of  the  different 
failure  criteria  as  they  pertain  to  anisotropic  materials. 

He  contrasts  the  use  of  the  maximum  stress  and  strain  crite- 
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ria  with  anisotropic  versus  isotropic  materials.  Wu 
(21:361)  noted  that  replacing  the  stress  invariants  with 
deviatoric  invariants  is  not  possible  since  hydrostatic  and 
deviatoric  yielding  has  not  yet  been  proven  to  be  indepen¬ 
dent.  Also,  as  shown  by  Jones  above,  material  parameters 
must  be  in  each  of  the  principal  directions  and  invariant  to 
transformations . 

A  more  complex  strength  criterion  was  proposed  by  Tsai 
and  Wu  to  take  into  account  stress  interactions  and  to  be 
invariant  in  material  transformations  (18:58-80).  This 
technique  uses  two  strength  tensors  in  a  scalar  function  to 
represent  a  failure  surface.  This  function 

/{^t)mFlal*Ft,ato r,-l  (2) 

contains  twenty-seven  terms  in  its  expanded  form  with  i,j,k 
=  1,2,... ,6.  The  coefficients  of  the  resulting  stress  poly¬ 
nomials  are  functions  of  the  material  strength  characteris¬ 
tics,  usually  determined  experimentally  (17:280-287).  The 
Tsai-Wu  criteria  has  favorable  properties  such  as  built  in 
invariance,  ease  of  transformation,  and  the  ability  to  show 
a  variety  of  stress  interactions  common  in  a  complex  field. 
But,  this  criteria  does  require  a  greater  knowledge  of  mate¬ 
rial  properties  and  can  become  somewhat  unwieldy. 


The  final  strength  criteria  to  be  considered  was  first  j 

< 

proposed  as  a  yield  criteria  for  anisotropic  materials  by  ] 

Hill  (6)  and  related  to  composite  strengths  by  Tsai  (16). 

This  Tsai-Hill  criteria  combines  the  apparent  simplicity  of 
maximum  strength  with  some  of  the  favorable  properties  of 
Tsai-Wu.  Principal  stresses  and  shear  stresses  are  combined 
in  the  following  form: 

I 

t 

1 

CQ  +  H)a2l  +  (F*H)al  +  (F  +  Q)a\-2H  a  ,a  2-2G<7  xa  3-2Fo  2a  3 

+  2Z.t|j  +  2A/t?3  +  2A/t?2-  1  (3) 

The  coefficients  F, G, H, L, M,and  N  are  functions  of  the  mate¬ 
rial  strengths.  These  properties  are  measured  experimental¬ 
ly  with  established  tension,  compression,  and  shear  tests. 

When  the  summation  of  the  combined  strength/stress  terms 
exceeds  unity  failure  is  predicted  to  occur.  According  to 
Jones  (8:79-80),  this  theory  works  well  with  E  glass/epoxy 
composites.  The  Tsai-Hill  criteria  also  reduces  correctly 
to  the  well  established  maxmum  octahedral  or  Von  Mises  shear 
stress  theory.  Tsai-Hill  can  be  regarded  as  an  accurate 
measure  of  composite  strength  for  the  load  cases  and  materi¬ 
al  orientations  being  investigated. 

Along  with  having  a  strength  criteria,  one  must  have  a 
way  of  analyzing  the  problem.  One  method  of  solving  complex 
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structural  problems  is  with  finite  elements.  "The  finite 
element  method  is  a  numerical  procedure  for  solving  a  con¬ 
tinuum  mechanics  problem  with  an  accuracy  acceptable  to  en¬ 
gineers"  (4:1).  With  the  advent  of  composite  materials  came 
the  application  of  finite  element  methods  to  these 
anisotropic  materials.  While  this  field  continues  to  grow 
as  the  "real  world"  is  modeled  more  closely,  significant 
achievements  have  already  occurred.  Some  of  these  achieve¬ 
ments  are  presented  here. 

Finite  element  modeling  of  composite  material  structures 
is  getting  more  sophisticated  as  people  try  to  use  this 
powerful  analysis  tool  in  more  applications.  Chamis  (2),  at 
the  NASA  Lewis  Research  Center  in  Cleveland,  Ohio,  reported 
on  their  current  work  using  finite  element  techniques  to 
analyze  composite  structures.  The  efforts  he  cites  include 
integrated  composite  mechanics  analysis,  simplified  compos¬ 
ite  mechanics  for  strength,  and  a  3-D  finite  element  for 
micromechanics  modeling  of  composites.  This  work  is  related 
closely  with  damage  and  damage  prediction.  Chamis  (2:5) 
refers  to  a  3-D  finite  element  model  "readily  adaptable  to 
investigating  interfacial  disbonds,  fiber  breaks,  voids  and 
even  matrix  crack  effects  on  the  microstress  distribution." 
This  "superelement"  is  not  yet  adaptable  to  many  finite 
element  routines,  but  serves  as  an  example  of  the  level  of 


sophistication  obtainable  with  finite  elements. 

Other  ways  of  using  the  finite  element  method  are  shown 
by  Ross  and  associates  (11).  They  use  SAPIV,  a  general 
purpose  code,  and  a  proprietary  code  named  ABAQUS  to  model 
interlaminar  delamination.  Anisotropy  is  accomplished  by 
varying  material  properties  along  the  different  axes.  ABAQUS 
has  an  interface  element  that  is  placed  between  brick  ele¬ 
ments  to  show  tnter laminar  behavior.  Their  results  indicate 
only  a  finite  progression  of  damage  since  they  are  modeling 
an  impact  event.  Delamination  is  the  failure  mode  they  are 
interested  in  studying,  Only  shear  stresses  at  discrete 
points  are  shown  as  the  authors  believe  that  shear  is  the 
major  cause  of  de laminat ion .  Humphreys  ( 7 )  uses  SAPIV  to 
develop  a  procedure  to  show  damage  accumulation  in  a  compos¬ 
ite  plate  during  a  low  velocity  impact.  He  calculates  dam¬ 
age  by  comparing  various  combinations  of  stress  and 
strength.  When  damage  is  said  to  occur  in  his  model  the 
failed  element  is  removed  from  the  mesh.  Humphreys  shows 
the  global  effects  of  this  impact  occurrence.  These  effects 
resemble  plate  theory. 

Damage  progression  work  done  by  Wi^t  (20)  used  a  two 
dimensional  finite  element  model  of  a  specimen  in  pure  ten¬ 
sion.  Failure  around  a  crack  was  predicted  using  two  dif¬ 
ferent  techniques.  The  first  method  used  a  purely  elastic 


technique  to  show  the  stress  field  after  damage  defined  by 
the  Tsai-Hill  criteria.  With  the  second  technique,  Witt 
would  completely  remove  the  stiffness  of  a  failed  element 
upon  damage.  While  neither  method  accurately  predicted  dam¬ 
age  growth,  they  did  provide  bounds  for  the  actual  phe¬ 
nomenon  (20:67)  . 

Shivakumar  and  Elber  (12)  used  the  finite  element  code, 
GAMNAS,  to  model  an  axisymmetric  plate  under  a  low  velocity 
impact.  By  using  the  crack  closure  technique,  they  sought 
to  identify  two  modes  of  delamination,  peeling  and  shear 
sliding.  They  only  used  cases  equivalent  to  static  deflec¬ 
tions  and  found  that  shear  sliding  was  the  primary  mode  of 
delamination  growth.  The  stress  fields  resulting  from  these 
cases  are  not  presented. 

Another  use  of  finite  elements  to  show  composite  damage 
is  by  Chamis  and  Williams  (3)  to  determine  interply  layer 
degradation  effects.  Applying  a  central  load  to  a  simply 
supported  beam,  they  looked  at  the  effects  that  a  reduced 
interply  (matrix  material)  stiffness  would  have  on  the  beam 
characteristics.  They  found  that  in  reducing  the  interply 
modulus  from  one  million  psl  to  one  thousand  psl,  no  notice¬ 
able  structural  changes  occur  until  around  the  ten  thousand 
ps i  level.  Although  not  showing  damage  progression  through 
the  beam,  they  do  show  some  overall  characteristics  that 


could  describe  a  beam  damaged  In  such  a  way  as  to  have 
reduced  stiffness  properties. 

Murthy  and  Chamis  (9)  used  a  NASTRAN  finite  element 
model  of  notched  and  smooth  fiber  composite  Charpy  specimens 
to  determine  failure  modes.  They  ran  both  static  and  dynam¬ 
ic  cases.  Their  results  showed  that  static  and  dynamic 
stress  contours  are  almost  identical.  The  quasi-static  rep¬ 
resentation  of  the  dynamic  process  was  the  most  conserva¬ 
tive.  While  not  concentrating  on  damage  progression,  they 
found  that  for  a  notched  specimen  either  normal,  longitudi¬ 
nal,  or  shear  stress  could  cause  or  initiate  a  failure. 

They  also  observed  that  shear  stresses  dominated  the  tran¬ 
sient  analysis  of  the  notched  specimen.  For  smooth  speci¬ 
mens,  combined  stresses  near  the  load  seemed  to  be  the  cause 
of  failure  with  no  one  stress  tending  to  dominate. 

Along  with  the  flexibility  of  finite  element  techniques 
come  certain  limitations.  Notably,  approximations  in  deter¬ 
mining  displacements  leads  to  further  approximations  in 
stress  calculations.  More  exact  methods  to  determine  the 
stress  field  using  numerical  techniques  or  integral  trans¬ 
forms  do  exist.  The  latter  technique  allows  one  to  solve 
the  mixed  boundary  value  problem  of  an  indenter  acting  on 
the  beam  surface.  These  closed  form  solutions  do  provide 
valuable  insight  with  which  to  compare  the  finite  element 
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analyses.  However,  they  can  become  rather  cumbersome  and 
can  not  be  readily  extended  to  include  damage  mechanisms. 

Shivakumar,  Elber,  and  Illg  (13)  use  numerical  methods 
to  solve  for  progressive  damage  in  thin  circular  laminates 
under  static  loads.  The  minimum  total  potential  energy 
method  combined  with  Von  Kantian  strain-displacement  equa¬ 
tions  gave  them  the  means  to  calculate  loads  as  well  as  ply 
stresses.  They  used  the  Tsai-Wu  criteria  to  determine  when 
failure  occurs  and  other  combinations  of  stresses  and 
strengths  to  evaluate  the  failure  modes.  Progression  is 
modeled  by  eliminating  the  stiffness  in  the  area  of  the 
calculated  failure  mode.  Their  analysis  concluded  that 
failure  started  at  the  bottom-most  ply.  In  terms  of  damage 
tolerances,  splitting  thresholds  were  lower  for  larger 
plates  while  "first-fiber-failure"  thresholds  were  lower  for 
smaller  plates. 

Another  way  of  solving  the  complex  mixed  boundary  con¬ 
tinuum  mechanics  equations  is  with  integral  transforms. 

Copp  and  Keer  (5)  analyzed  the  indentor  problem  with  simply 
supported  composite  beams  of  differing  fiber  orientations. 

An  extensive  parametric  study  of  contact  length,  layer 
thickness,  beam  length,  lamina  orientation  and  support  con¬ 
ditions  was  conducted.  The  resulting  stress  contours  de¬ 
picted  both  the  local  effects  of  the  indentor  and  the  sup- 


ports  as  well  as  the  global  beam  response.  While  not  using 
any  stiffness  reduction  scheme  for  damaged  material,  Copp 
and  Keer  did  use  the  Tsai-Hill  criteria  to  indicate  damage 
patterns  in  a  purely  elastic  material  over  various  load 
conditions.  The  integral  transform  techniques  used  did  not 
allow  for  a  progressive  increase  in  loading,  but  rather 
controlled  loading  by  specifying  the  contact  length  and  cal¬ 
culating  the  corresponding  indenter  loads.  Even  with  this 
disparity,  the  stress  contours  generated  by  Copp  and  Keer 
will  be  a  valid  comparison  for  any  stress  field  contours 
generated  by  this  thesis. 

Testing  provides  another  method  of  validation,  as  well 
as  a  means  to  explore  damage  progression.  Many  people  have 
tried  to  develop  empirical  damage  progression  models  based 
on  experimental  observations.  A  problem  with  using  test 
observations  is  that  damage  usually  occurs  quickly.  Even  if 
the  damage  progression  can  be  followed,  most  Instrumentation 
cannot  sort  out  the  complex  stress  field  closely  enough  to 
determine  which  stress  or  stresses  caused  the  failure.  For 
instance,  Whitney  and  Browning  (19)  point  out  some  apparent 
problems  with  the  well  established  short  beam  shear  test. 
They  claim  that  shear  is  not  the  most  dominant  of  the 
stresses  in  the  initial  stages  of  the  test.  The  first  mode 
of  failure  is  under  the  nose  of  the  load  applicator.  After 


initial  matrix  crushing  occurs,  shear  or  other  stress  combi¬ 
nations  then  complete  the  failure  process.  Their  results 
indicate  that  the  short  beam  shear  test  does  not  give  the 
expected  minimum  shear  strength  of  the  composite  material  in 
quest i on . 

Clearly,  modeling  the  phenomena  of  damage  initiation 
requires  a  failure  criteria.  Finite  elements  will  be  the 
method  of  calculating  stresses  and  applying  the  criteria. 
Finite  elements  will  also  be  the  means  of  incorporating 
stiffness  reduction  techniques.  The  stress  fields  will  be 
compared  to  those  of  Copp  and  Keer  (5)  for  accuracy.  Final¬ 
ly,  the  model  will  be  scrutinized  to  determine  how  closely 
it  matches  experimental  data. 


s 


'll 


This  thesis  Investigates  damage  initiation  in  a  beam 
under  a  contact  load.  Damage  initiation  is  defined  with 
respect  to  exceeding  a  damage  failure  criteria.  Without 
actually  describing  the  physics  of  the  damage  mechanism  such 
as  material  separation,  fracture,  or  plastic  deformation  one 
can  get  a  good  idea  of  where  and  how  the  material  would  fail 
by  looking  at  the  stress  fields  during  the  "failure"  pro¬ 
cess.  The  phenomenon  of  damage  initiation  as  it  pertains  to 
a  contact  load  on  a  beam  can  be  shown  with  a  typical  three 
point  bend  test  (Fig.  1).  The  short  beam  shear  test  de¬ 
scribed  by  ASTM  Standard  D2344  is  one  such  test.  With  this 
test  one  derives  apparent  shear  strength  properties  by  load¬ 
ing  the  specimen  to  failure  and  calculating  the  shear  by 


_  Q.75P, 

bd. 


(4) 


where  S  is  the  shear  strength,  ^  is  the  breaking  load  or 
maximum  load  on  the  load  indicating  mechanism  for  the  test, 
b  is  the  width  of  the  specimen,  and  d  is  the  thickness  of 
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the  specimen  (1:56-57).  This  test  is  relatively  simple. 


tempting  one  to  use  it  exclusively  in  establishing  design 


parameters  for  the  shear  strength  of  composite  materials 


But  as  seen  by  Copp  and  Keer  (5),  the  resulting  stress 


field  from  such  a  contact  problem  is  quite  complex,  thus 


failure  is  a  consequence  of  more  than  just  pure  shear.  A 


need  exists  to  determine  the  complete  state  of  stress  and  to 


translate  that  stress  field  into  an  estimation  of  damage 


Using  finite  elements  and  a  failure  criteria,  this  thesis 


analyzes  the  short  beam  shear  test  and  the  damage  initiation 


characteristics  of  different  fiber  orientations  in  parallel 


fiber  composites. 


Mbdel 


The  short  beam  shear  test  is  described  in  detail  by  ASTM 


Standard  D2344.  A  beam  is  supported  at  two  ends  with  the 


load  applied  in  the  center  (Fig.  1).  The  load  is  slowly 


increased  until  the  beam  specimen  breaks.  The  breaking  load 


is  defined  as  the  maximum  load  displayed  by  the  load  indi¬ 


cating  mechanism  used  for  the  test.  This  failure  load  is 


then  put  into  Eq  (4)  to  calculate  the  apparent  shear 


strength  of  the  specimen.  The  span  length- to -depth  ratio 


limits  for  test  specimens  is  five  (5)  for  materials  rein- 


forced  with  fibers  of  modulus  less  than  14.5  million  psl  and 
four  (4)  for  fibers  of  higher  modulus.  The  standard  also 
specifies  the  sizes  of  the  supports  and  the  load  nose.  This 
study  is  limited  to  the  basic  test  setup  and  specimen  ratios 
as  noted  in  the  standard.  The  loading  considered  is  only  a 
point  or  knife-edge  load. 

The  material  for  this  analysis  was  a  graphite/epoxy  type 
composite  with  a  Young's  modulus  fiber  to  resin  ratio  of 
about  four  (4).  The  material  properties  used  appear  in 
Table  1.  The  ratio  between  and  E ^  is  not  quite  as 

large  as  graphite/epoxy.  The  fourteen  to  one  ratio  usually 
associated  with  graphlte/epoxy  type  materials  could  not  be 
successfully  used  with  NOSAPM.  The  ratio  displayed  in  Table 
1  represents  the  largest  difference  in  moduli  completely 
checked  out  with  NOSAPM.  The  fibers  are  assumed  to  be  par¬ 
allel  with  material  properties  taken  as  the  average  over  a 
section.  Two  orientations  of  these  fibers  (Fig.  2)  are 
evaluated  to  determine  their  effects  on  damage  initiation. 
The  axis  along  the  fiber  is  represented  by  properties  denot¬ 
ed  by  ”11"  in  Table  1.  The  material  axes  perpendicular  to 
the  fiber  are  properties  "22"  and  "33".  One  would 


expect  that  the  ninety  (90)  degree  orientation  would  act 
most  like  isotropic  materials  because  of  the  in-plane  symme¬ 
try. 


Finite  Element  Model 

The  essential  parts  of  a  finite  element  analysis  include 
the  code  and  the  modeling.  This  research  employed  the  NOS- 
APM  code,  a  modified  version  of  NONSAP.  The  Aeropropulsion 
Laboratory  of  the  Air  Force  Wright  Aeronautical  Laboratories 
( AFWAL )  developed  NOSAPM  to  model  bird  strikes  on  turbine 
engine  blades. (15)  NOSAPM  is  a  multipurpose  finite  element 
code  capable  of  nonlinear  and  linear  static  and  dynamic 
analysis.  NOSAPM  uses  a  displacemental  formulation 
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where  the  stiffness  matrix  K  is 


[K}~  Vf  s[m]T E[m^B[m]cLV{ 
mJv 


with  B  dependent  on  the  shape  function  and  E  the  material 
properties  of  the  element.  The  vectors  U  and  R  are  the 


displacements  and  loads,  respectively.  A  one  dimensional 
truss  and  a  three  dimensional  quadrilateral  element  are 
available  in  NOSAPM.  The  element  properties  range  from  lin¬ 
ear  elastic  to  elastic-plastic.  A  linear  orthotropic  model 
exists  which  is  appropriate  for  representing  composite  mate¬ 
rials.  The  stress -stra in  relationship  for  this  orthotropic 
model  is 


where  C,  the  compliance  or  inverse  stiffness  matrix  is 


1  /£«  -vatt/Eb  -v./F, 

-vbJE%  l/Eb  -vbe/Ee 

l/E, 

0  0  0 

0  0  0 

0  0  0 

Note,  the  Young's  moduli,  shear  moduli,  and  Poisson's  ratios 
are  with  respect  to  the  material  axes.  Thus,  based  on  the 
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material  model  chosen,  NOSAPM  will  determine  the  displaceme¬ 
nts  first  and  then  calculate  the  stresses  based  on  these 
displacements.  In  addition  to  its  material  and  elemental 
capabilities,  NOSAPM  reformulates  the  stiffness  equations  at 
the  user’s  request  during  a  nonlinear  analysis  and  performs 
the  necessary  equilibrium  calculations  to  account  for  large 
displacements . 

As  mentioned  earlier,  the  short  beam  shear  test  was  the 
configuration  used  in  this  study  of  damage  Initiation.  In 
modeling  the  beam  for  the  short  beam  shear  test,  the  three 
dimensional  quadrilateral  element  was  used.  Each  of  these 
"brick”  elements  contained  eight  nodes,  one  for  each  corner 
of  the  element.  The  coordinate  system  chosen  uses  "x"  for 
the  thickness,  "y"  for  the  length,  and  "z"  for  the  depth  of 
the  beam,  as  well  as  for  the  direction  of  the  load  (Fig.  1). 
With  this  study  only  two  dimensions,  y  and  z,  were  left 
unconstrained  thus  creating  a  plane  strain  case.  With  a 
displacemental  solution,  this  restriction  allows  movement  in 
only  two  directions  reducing  the  three  dimensional  element 
to  two  dimensions.  Note  that  this  plane  strain  configura¬ 
tion  will  result  in  higher  stresses  than  the  plane  stress 
cases  of  Copp  and  Keer  (5)  used  for  validation.  The  basic 
shapes  of  the  stress  contours  should  be  sufficiently  close 
enough  for  the  comparisons. 


The  final  mesh  geometry  representing  the  beam  is  shown 
in  Figure  3.  As  shown,  the  model  consists  of  fourteen  (14) 
elements  through  the  thickness  with  eleven  (11)  elements 
along  the  beam  axis  for  a  total  of  fifty-four  (54)  elements 
The  model  is  only  one  element  thick.  Since  eight  noded 
brick  elements  are  employed,  the  total  number  of  nodes  is 
360.  A  convergence  study  determined  that  the  element  mesh 
should  be  refined  directly  under  the  load  in  order  to  accu¬ 
rately  model  the  strong  stress  gradients  in  this  area.  The 
model  takes  advantage  of  symmetry  along  the  centerline  of 
the  beam,  under  the  load.  Symmetry  is  required  because  of 
model  size  limitations  of  the  code  and  a  desire  to  conserve 
computer  resources.  The  support  is  modeled  as  a  constraint 
in  the  "z"  direction,  while  the  symmetry  nodes  along  the 
edge  are  constrained  from  moving  in  the  "yM  direction.  All 
nodes  are  constrained  in  the  "x"  direction. 


Damage  Cxltexla 

With  the  model  established,  the  next  step  is  to  incorpo 
rate  the  damage  criteria  into  the  computer  code.  Several 
different  criteria  were  considered  that  have  been  demon¬ 


strated  to  model  failure  or  damage  in  a  composite  material 
under  simple  loads.  These  included  maximum  stress. 


Tsai-Hill,  and  Tsai-Wu.  With  the  beam  test  configuration, 
the  state  of  stress  will  be  three  dimensional  and  interdepe¬ 
ndent.  The  failure  criteria  must  consider  the  influences  of 
all  the  stress  components.  Thus,  either  the  Tsai-Hill  or 
the  Tsai-Wu  appear  to  be  the  most  attractive  candidates. 

The  Tsai-Hill  theory  was  chosen  for  its  simplicity  and  ease 
of  adaptation  to  the  beam  test  problem.  Also,  with  the 
lamina  orientation  in  only  two  directions  of  zero  (0)  and 
ninety  (90)  degrees,  the  Tsai-Wu  theory  would  reduce  to 
Tsai-Hill  theory  for  this  orientation  and  geometry.  The 
Tsai-Hill  equations  are  found  in  Jones  (8:76-79)  as  follows: 


(G*  H)af+(F  +  H)o\+(F  +  G)a\  -  2H  a  xa  2-  2Go  ,a  3-2F  a  2o  3 

+  2Lx&  +  2M  X*3  +  2N  x*2-  1  (9) 


The  strength  parameters  are  F,  G,  H,  L,  M,  and  N  which 
relate  to  material  strength  properties  X,  Y,  Z,  and  S. 
These  relationships  are: 


with  S  being  the  shear  strength  in  the  respective  direc¬ 
tions. 


G  +  tf  - 


on 


Note,  X,  Y,  and  Z  are  the  tensile  or  compressive  strengths 
along  the  respective  material  axes.  The  remaining  three 
terms  for  the  cross  product  relationships  are: 


aHmh*h~h 


2G  -  -^r+  - - — 

X2  Z2  Y2 


2 F  -  -^+- - — 

Y2  Z2  X2 


(12) 


For  the  purpose  of  this  analysis  the  tensile  and  compressive 
strengths  of  each  axis  are  assumed  equal. 

Along  with  a  correct  determination  of  the  stress  fields 
the  damage  criteria  is  critical  in  the  study  of  damage  ini- 


tlation.  Damage  Initiation  begins  at  the  point  when  the 
damage  criteria  is  exceeded.  Subsequent  work  is  not  Intend¬ 
ed  to  analyze  the  stress  and  displacement  fields  in  the 
immediate  vicinity  of  the  damage.  Instead,  the  analyses 
reflect  the  containment  of  the  damaged  zone  by  the  surround- 
' ng  elastic  structure.  The  progression  of  the  damage  zone 
into  this  elastic  field  is  one  of  the  principle  interests  of 
this  research. 

The  Tsai-Hlll  failure  criteria  is  incorporated  into  NOS- 
APM  through  the  subroutine  DAMAGE.  A  discussion  and  docu¬ 
mentation  of  this  subroutine  is  given  in  Appendix  2. 

Briefly,  the  material  strength  properties  are  input  along 
with  the  other  material  characteristics.  A  control  card 
parameter  determines  which  criteria  will  be  used.  Current¬ 
ly,  DAMAGE  includes  only  the  Tsai-Hill  criteria.  It  should 
be  noted  that  any  desired  criteria  could  be  easily  added. 

As  NOSAPM  calculates  the  stresses  for  an  element  with  each 
step,  DAMAGE  checks  the  element's  stress  field  against  the 
criteria.  If  the  stress  field  exceeds  the  criteria,  the 
element  is  considered  damaged  and  is  so  noted  with  another 
control  parameter  in  NOSAPM.  Elements  are  checked  against 
the  criteria  individually  and  designated  as  damaged  accord¬ 
ingly.  A  damaged  element  can  then  be  modified  as  appropri¬ 
ate  using  reduced  stiffnesses  or  some  other  mechanism. 


Ift- 
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As  a  material  fails,  its  resistance  to  deformation  and 


its  ability  to  carry  loads  decreases.  One  way  this  can  be 


characterized  is  by  a  reduction  in  the  material’s  elastic 


properties.  This  reduction  is  intended  to  model  the  materi 


al  cracking  and  crushing  that  occurs  locally  but  precedes 


the  catastrophic  failure  of  de laminat i on .  With  the  finite 


element  approach  taken,  an  element  was  considered  failed 


when  the  stresses  within  the  element  combined  to  exceed  the 


Tsai-Hill  criteria  described  earlier.  While  this  thesis 


concentrated  on  the  initiation  of  damage,  an  insight  into 


the  progression  of  such  damage  is  possible.  By  tracking 


element  failures  and  making  the  appropriate  stiffness  modi¬ 


fications  one  can  discern  a  damage  growth  pattern.  This 


pattern  qualitatively  describes  the  damage  and  can  be  ex¬ 


plained  by  examining  the  surrounding  stress  fields.  Al¬ 


though  the  damaged  elements  do  not  accurately  reflect  the 


singular  stresses  present  due  to  cracks  or  plastic  deforma¬ 


tion,  the  results  do  coincide  with  experimentally  observed 


damage . 


Stiffness  reduction  is  accomplished  by  changing  the  ma¬ 


terial  properties  of  the  element.  With  the  input  file  sev¬ 


eral  sets  of  material  properties  can  be  specified  to  de¬ 


scribe  the  successive  decreases  in  stiffness.  The  first  set 


defines  the  material  properties  In  the  undamaged  state  The 


next  set  of  properties  Is  a  reduction  of  the  first  set  by 


some  predetermined  percentage.  The  following  sets  can  be 


further  reduced  as  desired.  Since  the  stiffness  is  directly 


determined  from  the  material  properties,  it  decreases  as  the 


properties  decrease.  A  control  variable  within  NOSAPM  calls 


a  new  set  of  properties  after  damage  occurs  and  with  each 


subsequent  stiffness  reformation.  The  result  is  an  element 


triggered  by  exceeding  some  failure  criteria  reducing  in 


stiffness.  As  the  loading  continues,  the  element  stiffness 


can  continue  to  be  reduced  through  the  use  of  additional 


sets  of  material  properties. 


Stiffness  reduction  was  implemented  two  different  ways 


The  first  was  a  reduction  which  occurred  with  each  step  in 


the  program.  The  stiffness  was  gradually  reduced  to  model 


an  increase  in  damage  to  elements  previously  failed.  Thus, 


at  any  step  the  element(s)  which  failed  first  would  be  soft¬ 


er  than  elements  just  failing.  The  second  method  of  stiff¬ 


ness  reduction  employed  was  to  impose  the  entire  reduction 


at  once.  Any  increase  in  program  steps  would  have  no  addi¬ 


tional  effect  on  failed  elements.  This  method  was  employed 


to  model  complete  crushing  and  cracking  of  an  element  con¬ 


strained  only  by  the  surrounding  elastic  field. 
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The  beam  bending  problem  analyzed  represents  a  stress 
field  with  a  combination  of  local  and  global  stresses.  The 
local  stresses  are  a  result  of  the  point  load  and  the  sup¬ 
ports.  The  global  stresses  are  the  beam  bending  stresses 
contained  in  any  beam  problem.  The  combinations  of  these 
stresses  determine  the  failures.  Three  stresses,  normal, 
bending,  and  shear,  are  important  in  this  two  dimensional 
representation  of  the  short  beam  shear  test.  These  stresses 
shown  in  Figure  4,  will  then  be  combined  through  the 
Tsai-Hill  relationships  to  determine  damage/ fa i lure . 

NOSAPM,  like  most  finite  element  programs,  calculates 
the  stresses  in  each  element  by  first  determining  dis¬ 
placements.  The  displacements  then  become  strains  and  Eq 
(7)  calculates  the  stresses  using  these  strains.  The 
stresses  as  well  as  displacements  are  then  tabulated  by  the 
program  with  the  output.  A  large  mesh  with  many  elements 
produces  a  large  table  of  stresses  for  each  step  in  the 
analysis.  These  tables  are  useful,  but  cumbersome. 

Post-processing  is  usually  required  to  transform  the  raw 
data  into  a  more  usable  form.  For  this  analysis  stress 
tuntours  of  the  stresses  calculated  at  Gauss  points  in  each 
element  were  plotted  by  a  program  called  SURFER,  a  graphics 
program  designed  to  produce  contours  and  surface  plots. 
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NOSAPM  required  some  modifications  to  produce  files  compati¬ 
ble  with  SURFER  so  that  it  could  be  utilized  as  a  post-pro¬ 
cessor.  Appendix  1  describes  these  NOSAPM  modifications  and 
the  resulting  output  files.  With  the  contour  plots,  the 
pages  of  stresses  are  reduced  to  single  pictures  for  rapid 
data  reduction. 

Each  stress  contour  is  important  in  the  determination  of 
damage  progression  through  the  short  beam  shear  model. 
Failure  mode  determinations  from  testing  usually  assume  the 
absence  of  some  stresses  or  the  dominance  of  other  stresses. 
For  example,  a  typical  tensile  strength  test  eliminates  all 
stresses  but  tension.  The  short  beam  shear  test  assumes 
that  the  transverse  shear  stress  dominates  and  results  in  a 
delamination  failure  defining  interlamina  shear  strength 
properties.  Finite  elements  will  allow  for  calculation  of 
each  stress  component.  Modifications  to  NOSAPM,  mentioned 
earlier,  permit  these  stresses  to  be  combined  by  the 
Tsai-Hill  criteria  to  define  failure.  Failure  can  then  be 
characterized,  and  the  dominant  stresses  identified.  Only 
then  can  a  failure  mode  be  deduced  and  the  damage  progres¬ 
sion  displayed.  Normal,  bending,  and  shear  stress  contours 
are  shown,  as  well  as  the  Tsai-Hill  stress  field. 
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Basalts  and  Else as si  an 

The  finite  element  program,  NOSAPM,  combined  with  the 
graphics  program,  SURFER,  produced  contour  plots  to  evaluate 
the  short  beam  shear  test.  Stress  fields  were  developed  to 
compare  with  other  more  exact  solutions.  Damage  zones  plot¬ 
ted  were  based  on  the  Tsai-Hill  failure  criteria.  The  con¬ 
tours  also  helped  to  evaluate  the  stiffness  reduction  tech¬ 
niques  as  a  method  of  determining  damage  initiation.  The 
results  below  attempt  to  satisfy  the  objectives  of  obtaining 
the  complex  stress  fields  of  this  mixed  boundary  value  prob¬ 
lem,  applying  a  failure  criteria  to  the  problem,  and  using 
stiffness  reduction  to  model  damage  initiation.  It  should 
be  noted  that  with  each  contour  plot  only  a  small  portion  of 
the  beam  is  presented  as  shown  in  Figure  5.  This  section  is 
under  the  load  and  along  the  line  of  symmetry  in  the  comput¬ 
er  model.  In  this  region,  the  local  stresses  are  most  sig¬ 
nificant.  Also,  this  region  contains  occurrences  of  ob¬ 
served  local  failures.  The  finite  element  mesh  has  been 
optimized  to  model  the  stress  gradients  in  this  region  and 
to  conserve  computer  resources . 

Shear  Test  Stress  Fields 

The  first  objective  of  this  study  was  to  determine 
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the  stress  fields  of  the  short  beam  shear  test.  High  stress 
gradients  were  expected  around  the  area  of  load  application. 
Thus,  the  greatest  refinement  of  the  element  mesh  should  be 
in  this  area.  A  convergence  study  demonstrated  this  need 
for  refinement  and  resulted  in  the  mesh  shown  in  Figure  3. 
The  mesh  emphasizes  the  outer  fibers  of  the  beam  as  well  as 
the  area  directly  below  the  applied  point  load.  Both  the 
ninety  degree  and  the  zero  degree  fiber  orientations  were 
used.  The  resulting  stress  contours  show  patterns  similar 
to  those  obtained  by  Copp  and  Keer  (5)  using  a  more  exact 
numerical  technique. 

The  normal  stresses  for  both  the  ninety  and  zero  degree 
cases  are  shown  in  Figures  6  and  7.  As  indicated  by  Copp 
and  Keer  (5:123,124),  these  stresses  are  dominated  by  the 
point  load  and  the  support  reaction.  Although  not  seen  with 
the  beam  section  shown,  these  normal  stresses  for  both  nine¬ 
ty  and  zero  degree  cases  disappeared  in  the  center  region  of 
the  beam  away  from  the  load  and  reaction.  Note  that  the 
normal  stress  contours  for  the  zero  degree  case  do  appear  to 
be  tighter  in  the  vicinity  of  the  load  point  than  the  ninety 
degree  case.  However,  the  deviation  is  small. 

Shear  values  for  the  same  cases  are  shown  by  Figures  8 


vary  from  the  load  point  to  the  center  line  as  expected  from 
beam  theory  and  observed  by  Copp  and  Keer  (5:126,127).  Un¬ 
like  the  normal  stresses,  the  shear  stresses  closely  agree 
for  the  two  cases.  These  shear  contours  were  the  most  dif¬ 
ficult  to  obtain.  Away  from  the  refined  mesh  area  these 
contours  wavered  and  became  ill-defined.  While  the  dis¬ 
placements  were  compatible  across  element  borders,  the  shear 
stresses  across  these  boundaries  were  not.  The  high  shear 
gradients  around  the  load  required  a  very  fine  mesh  to  ob¬ 
tain  stress  contours  as  detailed  as  Copp  and  Keer. 

Bending  stress  contours,  on  the  other  hand,  were  rela¬ 
tively  easy  to  generate.  The  ninety  and  zero  degree  bending 
cases  are  seen  in  Figures  10  and  11.  The  neutral  axis 
appears  along  the  centerline  of  the  beam  as  expected.  Al¬ 
though  not  shown  with  these  sections,  the  stresses  appear 
symmetric  about  a  neutral  axis  with  a  little  distortion 
around  the  point  load,  also  matching  the  patterns  found  by 
Copp  and  Keer  (5:129,130).  Beam  theory  clearly  predicts  the 
symmetry  about  a  neutral  axis.  The  distortion  around  the 
point  load  area  is  easy  to  understand  since  the  local  dis¬ 
placements  caused  by  the  loading  are  not  linear  with  respect 
to  the  rest  of  the  beam's  displacement. 


mesh,  the  stress  contour  patterns  seem  to  be  in  agreement 
with  the  work  of  Copp  and  Keer  (5).  A  similar  but  dynamic 
load  case  was  studied  by  Murthy  and  Chamis  (9)  resulting  in 
contour  patterns  close  to  those  shown  here.  It  now  seems 
reasonable  that  NOSAPH  along  with  the  mesh  used  can  deter¬ 
mine  the  complex  stress  fields  associated  with  the  short 
beam  shear  test  around  the  area  of  loading. 

Damage  Criteria  in  NOSAPM 

The  Tsai-Hill  criteria  was  used  as  the  means  of  deter¬ 
mining  elemental  failure.  The  subroutine  DAMAGE  (see  Ap¬ 
pendix  2)  compared  the  current  stress  fields  in  each  element 
with  the  material's  strengths.  Damage  occurs  when  the  com¬ 
binations  of  stresses  and  strength  ratios  are  greater  than 
one  (1).  The  resulting  plots  for  both  the  ninety  and  zero 
degree  cases  show  a  state  of  stress  past  initial  failure. 
Figures  12  and  13  show  stress  combinations  exceeding  the 
Tsai-Hill  criteria  with  values  greater  than  one.  The  over¬ 
load  case  may  be  somewhat  artificial  in  the  sense  that  mate¬ 
rial  degradation  would  already  have  occurred  and  the  failure 
pattern  would  probably  experience  some  change.  However,  the 
overload  cases  are  a  good  comparison  to  the  works  of  Copp 
and  Keer  (5:140,141). 
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The  ninety  degree  case  represents  material  that  is 
isotropic,  especially  in  the  plain  strain  configuration 


studied.  The  failure  pattern  shown  (Figure  12)  closely  re¬ 
sembles  the  bending  stress  contours  (Figure  8).  This  pat¬ 
tern  shows  some  localized  effects  from  the  load,  but  primar 
ily  the  failure  exhibited  is  beam  bending.  The  zero  degree 
case  is  representative  of  anisotropic  or  composite  materi¬ 
als.  Its  failure  pattern  (Figure  13)  most  closely  matches 
the  normal  and  shear  stress  cases  (Figures  9  and  10).  This 
indicates  that  the  predominant  failure  modes  are  normal  and 
shear  with  the  load  application  point  the  major  contributor 
The  absence  of  the  beam  bending  patterns  in  this  case 
suggests  that  normal  beam  theory  is  not  applicable.  The 
extremely  high  failure  values  near  the  load  point  seem  to 
reinforce  the  theory  by  Whitney  and  Browning  (19)  that  the 
short  beam  shear  test  possesses  a  more  complex  stress  field 
than  beam  theory  predicts.  The  next  step  in  this  damage 
initiation  study  was  to  model  damage  by  reducing  the  stiff¬ 
ness  in  the  failed  element(s). 

St iJLEnss?  Reduction  la  NOSAPM 

Stiffness  reduction  was  incorporated  into  NOSAPM  using 
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the  subroutine  DAMAGE  and  the  code's  ability  to  handle  a 
number  of  different  sets  of  material  constants.  Property 
sets  were  read  into  NOSAPM  with  the  moduli  reductions  corre¬ 


sponding  to  desired  stiffness  reductions.  For  example,  if  a 
twenty  percent  reduction  was  needed,  the  first  set  of  mate¬ 
rial  properties  would  be  the  originals  and  the  second  set 
would  be  values  reduced  by  twenty  percent.  The  Poisson's 
ratio  would  remain  the  same,  while  the  Young's  modulus  and 
shear  modulus  would  be  reduced  by  twenty  percent.  When  DAM¬ 
AGE  determined  a  "failed"  element,  NOSAPM  would  then  assign 
the  second  set  of  properties  to  that  element  while  the  other 
elements  retained  their  original  values. 

This  stiffness  reduction  technique  affects  the  entire 
element  once  failure  occurs.  First,  the  technique  was 
checked  with  a  simple  single  element  compression  case. 

Then,  the  technique  was  used  with  the  short  beam  shear  test 
model  with  fibers  oriented  in  both  the  zero  and  ninety  de¬ 
gree  configurations.  Two  questions  are  addressed  with  this 
study.  The  first  is  whether  the  technique  works  with  the 
finite  element  code  NOSAPM.  The  second  question  deals  with 
the  applicability  of  modeling  damage  initiation  by  reducing 
the  stiffness  of  the  entire  element. 

Single  Element  Compress  ion .  As  a  test  of  the  stiffness 
reduction  modifications  to  NOSAPM  a  single  element  under 
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simple  compression  was  investigated.  The  element  was  con¬ 


strained  on  three  sides  with  a  compressive  load  in  the  "z" 
direction  (Figure  14).  The  load  was  held  constant  while  the 
stiffness  was  reduced  five  percent  (5  %)  with  each  solution 
step.  Simple  mechanics  theory  provided  the  predictions  for 
d isplacements  and  stresses. 

To  begin  with,  the  displacements  for  the  single  element 
can  be  found  using  the  axial  deformation  equation: 


fi- 


PL 

AE 


(13) 


where  P  is  the  load,  L  is  the  original  length,  A  is  the 
original  area,  and  E  is  the  Modulus  of  Elasticity.  The  load 
remained  constant  during  the  entire  run.  The  length  and 
area  can  be  considered  constant  since  the  changes  caused  by 
the  load  considered  are  small  with  respect  to  the  original 
dimensions.  With  the  above  assumptions,  one  would  then  ex¬ 
pect  the  deflection  to  be  inversely  proportional  to  the 
modulus . 
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Another  parameter  to  check  with  the  program  is  the 


stress.  In  this  simple  configuration,  the  stress  can  be 
approximated  by  the  axial  stress  equation: 


As  mentioned  earlier,  the  load  is  constant  and  the  area  is 
assumed  unchanged.  The  stress  within  the  single  element  can 
then  be  expected  to  remain  constant  as  it  is  independent  of 
any  modulus  change. 

The  results  of  this  check  case  on  the  behavior  of  the 
stiffness  reduction  scheme  are  shown  in  Table  2.  The  calcu¬ 
lated  displacements  do  not  vary  from  theory  except  for 
slight  round-off  and  truncation  differences.  The  only 
stress  possible,  that  in  the  "z"  direction,  remained  con¬ 
stant  as  predicted.  These  results  give  confidence  that  the 
stiffness  reduction  scheme  does  not  interfere  with  the  dis- 
placemental  solutions  of  NOSAPM  and  correctly  computes 
stresses  based  on  the  most  current  stiffness  moduli  values. 

Zero  Degree  Fiber  Orientation .  The  effects  of  stiffness 
reduction  as  a  possible  model  of  damage  initiation  in  com¬ 
posites  was  studied  using  a  short  beam  shear  test  configura- 
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tlon  developed  earlier  with  fibers  in  the  zero  degree  orien¬ 
tation.  This  parametric  study  included  two  extremes.  The 
first  was  a  constant  load  with  a  varying  reduction.  The 
second  increased  the  load  while  keeping  the  reduction  fixed. 
These  two  cases  provided  bounds  for  this  reduction  tech¬ 
nique  . 


Single  Element  Check 
for 

Stiffness  Reduction 
Table  2 
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As  stated  earlier,  the  first  case  kept  a  constant  load 
with  a  varying  stiffness  reduction.  This  constant  load  of 
three  (3)  pounds  was  chosen  as  it  represented  the  first 
indication  of  damage  occurrence.  NOSAPM  showed  the  element 
directly  under  the  load  failing  at  the  integration  point 
nearest  to  the  load.  The  output  file,  TAPE30,  was  then  used 
to  track  any  changes  to  element  failure  with  stiffness  re¬ 
ductions.  TAPE30  and  other  input/output  files  are  explained 
further  in  Appendix  A.  Stiffness  reductions  produced  the 
threshold  of  damage  initiation  in  the  element.  The  progres¬ 
sion  of  damage  into  the  body  can  only  be  observed  by  noting 
the  increase  in  the  stress  fields  of  the  adjacent  elements 
as  shown  in  Figures  15,  16,  and  17.  From  100%  (Figure  15) 
to  about  50%  of  the  original  stiffness,  the  failed  element 
does  not  change  as  stresses  at  only  one  integration  point 
exceeded  the  Tsai-Hill  criteria.  The  damage  contour  shows 
the  dominance  of  the  normal  stress.  Between  50%  (Figure  16) 
and  about  5%  of  the  original  stiffness,  stresses  at  two 
integration  points  within  the  element  exceed  the  Tsai-Hill 
criteria.  Contour  lines  tighten  up  around  the  load  but  seem 
to  spread  out  more  through  the  beam.  A  stiffness  reduction 
smaller  than  5%  (Figure  17)  of  the  original  stiffness  trans¬ 
fers  the  load  through  the  element  in  such  a  way  that 
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Eq  (9),  Tsai-Hill,  0%  Stiffness  Reduction 
Constant  Load,  P  =  3  lb 
Fig.  15 
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no  failure  occurs  within  the  element.  In  other  words, 
stresses  within  the  reduced  stiffness  element  at  all  inte¬ 
gration  points  fall  below  the  failure  criteria  of  Tsai-Hill. 
The  intensity  of  the  applied  load  as  shown  in  the  Tsai-Hill 
plots  appears  to  be  spread  into  the  beam,  reducing  the 
stresses  in  the  originally  failed  element. 

The  normal  and  shear  stress  contours  remain  the  same 
shape  regardless  of  stiffness  reduction.  Some  differences 
can  be  seen  from  the  full  stiffness  case  (Figures  18  and 
19).  Stiffness  reduction  appears  to  affect  the  normal 
stresses  by  shifting  them  away  from  the  load  point  as  shown 
in  Figures  20  and  21.  The  maximum  normal  stress  is  lower 
with  the  larger  reduction,  but  higher  stresses  are  spread 
further  in  the  beam  (Figure  21).  The  same  is  true  for  the 
shear  stress.  With  the  fifty  percent  reduction  case  (Figure 
22),  the  shear  increases  around  the  load.  As  the  stiffness 
is  reduced  further  to  five  percent  of  its  original  stiffness 
(Figure  23),  the  shear  decreases  around  the  load  while  in¬ 
creasing  in  the  surrounding  elements.  The  increasing  shear 
seems  to  be  moving  towards  the  centerline  of  the  beam.  This 
behavior  does  correlate  with  observed  patterns  of  failure  in 
short  beam  shear  tests  (20).  Maximum  shear  is  expected 
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along  the  centerline,  but  the  observed  failure  path  seems  to 
originate  close  to  the  load  and  run  to  the  centerline.  The 
diagonal  path  taken  is  similar  to  that  shown  by  these  shear 
stress  contours. 

The  biggest  differences  caused  by  the  stiffness  reduc¬ 
tion  technique  occurs  in  the  bending  stresses.  Figure  24 
shows  the  contours  resulting  from  the  original  single  ele¬ 
ment  damage.  The  reduced  stiffness  cases  are  the  same  when 
looking  away  from  the  damaged  element.  However,  stiffness 
reduction  changes  the  bending  stresses  locally.  The  fifty 
percent  case  (Figure  25)  shows  the  compression  lines  moving 
toward  the  load  on  top  of  the  beam  while  moving  away  from 
the  load  through  the  beam's  thickness.  As  the  stiffness  is 
reduced  to  five  percent  (Figure  26),  the  point  of  maximum 
compression  shifts  away  from  the  surface  to  an  area  under 
the  damaged  element.  The  reduced  stiffness  element  seems  to 
relieve  some  compression  while  shifting  the  load  down.  This 
dramatic  change  in  bending  stress  patterns  is  not  too  sur¬ 
prising.  The  zero  degree  case  shown  has  the  highest  modulus 
along  the  axis  of  bending.  Reduction  causes  a  larger  dif¬ 
ference  in  modulus  along  this  axis  than  for  the  normal  or 
shear  cases.  This  big  modulus  difference  then  changes  the 
distribution  of  bending  stresses  as  shown.  Although  reduced 


the  same  percentage,  the  differences  in  the  other  material 


property  values  are  not  as  great.  So  normal  and  shear 


stresses  do  not  change  as  dramatically.  Since  the  failure 


is  dominated  by  normal  and  shear,  the  Tsai-Hill  plots  would 


not  change  dramatically  with  reduction. 


The  next  limiting  case  tried  with  a  zero  degree  fiber 


orientation  was  to  increase  the  load  while  holding  the  dam¬ 


aged  element  reduction  constant  at  five  percent  of  the  orig¬ 


inal  value.  The  load  was  a  ramp  increase  from  five  to  twen¬ 


ty  pounds  in  four  steps.  The  program  ran  for  only  two  steps 


and  stopped  with  a  numerical  instability.  By  the  end  of  the 


second  step,  displacements  in  the  failed  element  were  very 


large.  Also,  during  the  third  step  a  non-positive  definite 


stiffness  matrix  was  formed  causing  program  termination. 


The  bending  stress  (Figure  27)  is  similar  to  the  constant 


load  case  previously  discussed.  The  normal  (Figure  28)  and 


shear  (Figure  29)  stresses  show  maximums  shifting  from  the 


surface  into  the  beam.  The  normal  stress  shows  a  zone  of 


highest  compression  inside  the  beam.  Physically  this  seems 


unreasonable  as  an  area  of  tension  cannot  exist  directly 


under  a  compressive  load.  The  damaged  area  seems  to  possess 


an  unloading  action  resulting  from  the  element  stiffness 


reduction.  This  is  seen  as  the  maximum  compression  appears 
to  move  into  the  beam.  The  shear  also  goes  from  negative  to 
positive  with  the  high  load.  This  is  due  to  the  high  stress 
gradients  and  the  model  not  having  enough  detail  to  pick  up 
all  the  behavior  locally.  The  continuation  of  loading  with 
stiffness  reduction  activated  appears  to  result  in  the  load 
moving  through  the  beam.  However, the  reversal  in  expected 
stresses,  i.e.  tension  instead  of  compression,  causes  this 
stiffness  reduction  technique  to  be  suspect. 

Ninety  Degree  Fiber  Or ientation .  Along  with  the  zero 
degree  case,  a  ninety  degree  case  was  run  to  consider  the 
"isotropic"  behavior.  The  load  was  held  constant  while 
stiffnesses  were  reduced.  Unlike  the  zero  degree  case,  re¬ 
ductions  through  five  percent  of  the  original  stiffness 
caused  no  noticeable  changes  in  damage  or  stress  contours. 
The  high  modulus  axis  was  along  a  direction  not  studied. 

The  axes  of  interest  had  moduli  of  similar  low  values. 
Apparently,  even  a  high  percentage  of  reduction  did  not 
separate  the  material  properties  significantly  to  change 
stresses  or  displacements.  The  constant  load  was  increased 
to  levels  far  higher  than  the  initial  failure  levels  with  no 
apparent  changes  due  to  stiffness  reductions.  With  no 
changes  in  the  constant  load  case,  an  increasing  load  with 
constant  reduction  case  was  not  attempted. 


The  short  beam  shear  test  in  composite  materials  was 
studied  with  finite  elements  to  determine  a  means  of  model¬ 
ing  damage  initiation.  The  finite  element  routine  used  was 
NOSAPM,  a  modified  version  of  NONSAP.  The  complex  stress 
fields  associated  with  this  shear  test  were  generated  suffi¬ 
ciently  around  the  area  of  the  load  application.  A  subrou¬ 
tine  to  apply  a  damage  criteria  was  incorporated  into  NOSAPM 
to  determine  element  "failure"  in  the  complex  stress  field. 
Also,  a  stiffness  reduction  technique  was  applied  where  the 
entire  failed  element  was  softened  by  varying  degrees.  The 
results  of  this  study  indicate  the  following: 

(1)  The  convergence  study  revealed  that  normal  and  bending 
stresses  can  be  obtained  fairly  well,  while  shear  stress 
is  a  function  of  the  mesh  size. 

(2)  The  Tsai-Hill  failure  criteria  seems  to  be  readily 
adaptable  to  the  finite  element  method. 

(3)  Fibers  oriented  ninety  degrees  appear  to  be  unaffected 
by  the  stiffness  reduction  technique  especially  with  the 
low  modulus  properties  used. 


(4)  For  zero  degree  fiber  orientation,  reductions  up  to 
fifty  percent  result  in  no  drastic  changes  in  failure  or 
stress  contours. 

(5)  Reductions  between  fifty  and  ninety-five  percent  for 
zero  degree  fibers  show  more  damage  for  a  constant  load. 

(6)  Stiffness  reductions  above  95%  for  the  zero  degree  fiber 
case  showed  the  originally  failed  element  no  longer 
exceeded  the  failure  criteria,  while  the  load  was  car¬ 
ried  by  the  adjacent  elements. 

(7)  Increasing  the  load  with  ninety-five  percent  stiffness 
reduction  in  zero  degree  fiber  orientation  cases  does 
not  change  the  pattern  of  failures  and  stress  contours 
except  to  move  the  maximums  into  the  body  and  away  from 
the  surface  being  loaded. 

To  further  this  study  of  damage  initiation,  one  should 
determine  the  exact  failure  modes  occurring  and  then  reduce 
the  stiffness  in  only  that  part  of  the  element.  For  exam¬ 
ple,  if  the  mode  of  failure  was  bending,  only  the  modulus 
along  the  bending  plane  would  be  reduced.  Load  shifting 
techniques  should  also  be  studied  to  circumvent  program  ter¬ 
mination  due  to  non-positive  stiffness  matrix  formulations 
on  ill-conditioned  elements.  Finally,  by  putting  a  model  of 
the  indenter  into  the  analysis  one  could  avoid  the  severe 
point  load  gradients  produced  in  the  single  element  case. 


The  NOSAPM  Users  Manual  (15)  describes  in  detail  the  use 
of  this  finite  element  code.  Modifications  were  necessary 
to  incorporate  the  DAMAGE  subroutine  outlined  in  Appendix  B 
and  the  post-processing  graphics  program  SURFER  (14).  This 
section  deals  with  three  classes  of  mod i f icat ions :  (1)  the 

input  files  for  NOSAPM,  (2)  the  job  control  listing  (JCL) 
for  the  Wr ight-Patterson  AFB  (WPAFB)  CYBER  computer,  and  (3) 
the  output  data  files  for  SURFER.  Discussions  will  include 
the  steps  to  running  both  the  main  program  and  the  associat¬ 
ed  post-processing. 

Pr.<?P3t ;nq  Input  Files 

The  NOSAPM  input  data  cards  are  clearly  explained  in  the 
users  manual  (15:1-87).  The  instructions  given  should  be 
followed  exactly  except  as  indicated  in  Figure  30.  Note 
that  changes  and  specific  values  needed  to  run  the  stiffness 
reduction  cases  are  underlined.  Also,  any  cards  not  shown 
in  the  figure  can  be  assumed  to  be  the  same  as  described  in 
the  users  manual.  As  can  be  seen,  changing  the  code  to  get 


the  necessary  formatted  data  to  SURFER  restricted  NOSAPM  to 
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the  use  of  rectangular  matrices  of  three  dimensional  ele¬ 
ments  using  a  "2X2"  Gaussian  integration  order.  By  using 
another  post-processing  technique  in  plotting  stresses  (or 
strains)  one  could  eliminate  these  restrictions. 

Running  MQSAPM  at  wright-Patterson  &££ 

For  this  study,  NOSAPM  was  run  on  the  CYBER  computer  at 
WPAFB .  The  length  of  this  program  requires  that  It  be  run 
as  a  batch  job  using  control  cards  or  a  JCL.  An  example  of 
a  typical  JCL  for  this  study  is  shown  in  Figure  31.  The 
steps  used  for  this  study  are  as  follows.  After  initiating 
the  batch  job,  certain  direct  access  files  are  opened  for 
processing.  The  file  titled  "NOSRUN"  is  the  compiled  ver¬ 
sion  of  NOSAPM.  The  three  other  input  files,  TAPE91, 
TAPE31,  and  TAPE30,  are  used  within  NOSAPM  to  store  SURFER 
formatted  data.  TAPE91  contains  the  bending,  normal,  and 
shear  stresses  matched  to  their  corresponding  integration 
points  on  the  face  of  the  beam.  These  stresses  are  normal¬ 
ized  to  the  material's  compressive  strength  values.  TAPE 31 
lists  the  failure  criteria  values  at  each  integration  point 


on  the  front  face  of  the  beam.  TAPE30  shows  the  elements 


JWFAK , T2  00  0 . 

USER, P870075, DUPOQE1 . 
CHARGE, * . 

GET , TAPE5=FNM0SM . 
ATTACH, NOSRUN. 

DEFINE, TAPE91. 

DEFINE, TAPE30 . 

DEFINE, TAPE31. 

NOSRUN, PL=80000 . 
RETURN, TAPE30 . 

RETURN, TAPE31 . 

RETURN , TAPES  1 . 

RETURN, NOSRUN. 

#EOR 

#EOI 


Job  Control  Listing  for  NOSAPM 
Fig.  31 
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which  have  exceeded  the  failure  criteria  and  the  integration 


point(s)  in  the  element  where  the  criteria  was  exceeded.  It 


is  important  to  remember  that  in  using  this  JCL  one  should 


not  already  have  files  TAPE91,  TAPE31 ,  and  TAPE30  in  perma¬ 


nent  storage.  The  JCL  creates  these  files  before  each  run. 


These  output  files  (TAPE91,  TAPE31,  and  TAPE30)  must  then  be 


removed  or  renamed  prior  to  subsequent  runs. 


Using 


for  Contoui 


With  the  output  files  from  NOSAPM,  contour  plotting  of 


stresses  and  failure  criteria  are  possible.  The  graphics 


package  used  was  SURFER ,  a  personal  computer  resident  plot¬ 


ting  program.  The  NOSAPM  output  files  had  to  be  downloaded 


to  a  floppy  disk  using  a  communications  package  with  an 


acceptable  protocol.  The  CYBER  had  XMODEM  protocol  avail¬ 


able  for  this  transfer  of  files  from  the  mainframe  computer 


to  a  personal  computer.  Each  NOSAPM  generated  file  usually 


contained  more  than  one  time  step.  So,  the  entire  listing 


for  each  time  step  had  to  be  separated  when  placed  on  the 


floppy  disk.  SURFER  does  accept  the  ASCII  character  format 


resulting  from  this  file  transfer. 


Actual  contouring  of  the  data  is  best  described  by  the 


SURFER  users  manual  (14).  Note  that  NOSAPM  generates  the 


data  for  the  entire  beam.  Since  only  a  small  part  of  the 
beam  is  important  in  this  study,  much  of  the  NOSAPM  data  wa 
truncated  prior  to  contouring.  This  truncation  saved  pro¬ 
cessing  time  and  allowed  for  a  more  exact  contour.  The 
contours  are  easily  viewed  with  SURFER  and  displayed  by  the 
computer  screen,  a  printer,  or  a  plotter. 
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The  purpose  of  subroutine  DAMAGE  was  to  incorporate 


failure  criteria  into  the  program  NOSAPM.  It  is  called  dur¬ 


ing  the  stress  calculation  portion  of  the  subroutine  THDFE . 


After  the  stresses  are  calculated  for  each  element  at  the 


integration  points,  THDFE  calls  DAMAGE  to  compare  the 


stresses  against  a  predetermined  failure  criteria.  So,  DAM¬ 


AGE  combines  the  stress  fields  generated  at  each  integration 


point  for  comparisons  against  strength  criteria.  The  sub¬ 


routine  is  listed  in  Figures  32a  &  32b. 


Input  to  DAMAGE  is  through  common  statements.  The  con¬ 


trol  character  for  failure  and  stiffness  reduction  is 


IFAIL(N)  where  N  is  the  element  number.  The  yield  stress 


and  maximum  strain  are  denoted  by  SSFAIL  and  SNFAIL,  respec¬ 


tively.  The  control  character  IFCRIT  tells  the  routine 


which  failure  theory  to  use.  Currently  DAMAGE  only  includes 


the  Tsai-Hill  criteria.  Different  theories  such  as  maximum 


stress  and  maximum  strain  can  easily  be  added  to  DAMAGE. 


Other  input  to  the  subroutine  are  the  stresses  or  strains  of 


the  element.  These  are  calculated  at  each  integration  point 


of  the  element  with  another  routine  in  NOSAPM  and  enter 


through  the  common  statement  MTMD3D. 


For  the  Tsai-Hill  theory  listed,  the  code  follows  the 


pr-.  izr  PL-’  V".  V'  mji  VV^.V-V'^  L^IF  [W1  'JV  LTV  'J*VV  V*  \TKV*V«  V*1T»  iT*  Y*  IT*  UT" 


p' 


s. 
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equations  outlined  by  Jones  (8).  DAMAGE  computes  failure 
using  stresses  and  strengths  in  combinations  determined  by 
the  Tsai-Hill  relationships.  These  combinations  within  DAM¬ 
AGE  result  the  dimensionless  parameter  (FTSAI)  which  if 
greater  than  one  (1)  indicates  exceedance  of  the  criteria. 
For  the  Tsai-Hill  case,  ETSAI  is  the  actual  dimensionless 
value  of  the  stress  combinations.  FTSAI  is  the  check  value 
and  is  equal  to  ETSAI  minus  one  (1). 

Another  calculation  performed  by  DAMAGE  is  to  check  the 
stiffness  reduction  control  character  IFAIL(N).  At  the 
start  of  NOSAPM,  IFAIL(N)  is  set  to  zero  (0)  for  each  ele¬ 
ment.  Upon  failure  of  the  element,  IFAIL(N)  becomes  one 
(1).  IFAIL(N)  increments  by  one  (1)  with  each  succeeding 
step  of  the  program.  With  this  method,  the  element(s)  fail¬ 
ing  on  earlier  time  steps  would  have  a  higher  IFAIL(N)  value 
than  those  elements  just  failing.  In  NOSAPM  I FAIL ( N )  de¬ 
notes  the  material  property  set  number  for  the  element.  The 
material  property  set  dictates  the  state  of  elemental  fail¬ 
ure  as  mentioned  in  Appendix  A. 

The  output  of  this  subroutine  is  tailored  to  the  program 
SURFER.  ETSAI  is  stored  to  create  contour  plots  of  the 
Tsai-Hill  values.  The  raw  data  collected  here  is  combined 
in  the  main  program  with  the  position  coordinates  of  the 
respective  integration  point  to  complete  the  package  for 
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SURFER.  Another  piece  of  the  output  set  up  for  this  study 
is  the  listing  of  the  failed  elements  and  the  integration 
points  where  failure  occurred.  Thus  one  can  get  a  quick  and 
easily  plotted  picture  of  damage  as  it  progresses  not  only 
through  the  integration  points  of  the  element  but  also 
through  the  entire  model. 
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SUBROUTINE  DAMAGE 
C 
C 

C  SUBROUTINE  TO  CHECK  FOR  FAILURE  IN  AN  ORTHOTROPIC 

MATL 

C 

C  IFCRIT  =  0  NO  FAILURE  CRITERIA 
C  =1  MAXIMUM  STRESS  USED 

C  =2  MAXIMUM  STRAIN 

C  =3  TSAI-HILL  THEORY 

C 
C 

COMMOM/FAIL/IFAIL ( 300 ) , SSFAIL( 3, 3 ) , SNFAI L ( 3 , 3 ) , IFCRIT 
COMMON/MTMD3D/D (6,6 ) , STRESS ( 6 ) , STRAIN ( 6 ) , IPT,N 

C 

C 

I F ( IFCRIT .EQ . 1 )  GO  TO  100 
I F ( IFCRIT. EQ. 2)  GO  TO  200 
I F ( IFCRIT. EQ.0)  RETURN 
C 
C 


CL 

= 

1 . 0 

/ 

( SSFAI L ( 3 , 3 ) *  *  2 ) 

CM 

=S 

1.0 

/ 

(SSFAIL(2,3)**2) 

CN 

= 

1 . 0 

/ 

{ SSFAIL ( 1, 3 ) **2 ) 

CA 

= 

4  .0 

/ 

( SSFAI L ( 1 , 1 ) *  *  2 ) 

CB 

1.0 

/ 

( SSFAIL ( 2, 1 ) **2 ) 

CC 

= 

1.0 

/ 

( SSFAIL ( 3, 1 ) **2 ) 

CD 

s 

CA+CB- 

■CC 

CE 

= 

CA+CC- 

CB 

CF 

= 

CC+CB- 

CA 

FTSAI  =CA* ( STRESS ( 1 ) **2  )  +CB* ( STRESS ( 2 ) * * 2 )  + 

&  CC* ( STRESS ( 3 ) *  *  2 )  -CD*STRESS ( 1 ) *STRESS ( 2 )  - 

&  CE*STRESS ( 1 ) *STRESS ( 3 ) -CF*STRESS ( 2 ) *STRESS ( 3 ) + 
&  CL* (STRESS(6)**2)  +  CM* ( STRESS ( 4 ) **2 )  + 

&  CN* ( STRESS { 5 ) **2  )  -  1.0 

ETSAI  =  FTSAI  +  1.0 
IF ( IPT . EQ . 5 )  WR I TE (96,4000)  ETSAI 
I F ( I PT . EQ . 6 )  WRITE(96, 4000)  ETSAI 
IF(IPT.EQ.7)  WRITE(96, 4000 )  ETSAI 
I F ( I PT . EQ . 8 )  WRITE( 96, 4000 )  ETSAI 
4000  FORMAT ( Fl 2 . 5 ) 

IF ( FTSAI . GT .0.0)  GO  TO  300 
RETURN 


Fig.  32a:  Source  Listing  of  Subroutine  DAMAGE 
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c 

c 

100  RETURN 
C 
C 

200  RETURN 
C 
C 

300  IF( IFAIL(N) . EQ . 0 )  IFAIL(N)  =  1 
IF ( I PT . LT . 8 )  GO  TO  301 

IF(IFAIL(N) .GT.O)  IFAIL(N)  =  IFAIL(N)  +  1 

301  WRITE( 30, 1000)  N,IPT 

1000FORMAT(8X,16H  ELEMENT  FAILED  , I4,20H  INTEGRATION 
POINT  ,14) 

RETURN 

END 


P 

P 

P 


Fig.  32:  Source  Listing  of  Subroutine  DAMAGE 
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The  purpose  of  this  study  was  to  analyze  damage  initia¬ 
tion  in  the  short  beam  shear  test  for  composite  materials. 
This  thesis  proposes  a  means  of  predicting  damage  initiation 
in  a  composite  beam  using  finite  element  techniques  and  a 
continuum  mechanics  failure  criteria.  The  first  objective 
was  to  accurately  determine  the  stress  field.  The  second 
objective  was  to  employ  a  damage  criteria  based  on  the  ele¬ 
mental  stress  field.  The  final  objective  was  to  evaluate 
the  effectiveness  of  using  an  elemental  stiffness  reduction 
to  study  damage  initiation  and  progression.  The  finite  ele¬ 
ment  code  used  was  NOSAPM. 

The  complex  stress  fields  associated  with  this  problem 
were  accurately  obtained  with  the  exception  of  the  shear 
stresses.  The  high  shear  gradients  required  the  use  of  a 
very  fine  mesh  which  was  not  possible  with  the  computer 
resources  available.  The  Tsai-Hill  failure  criteria  was 
successfully  incorporated  into  NOSAPM.  Also,  an  elemental 
stiffness  routine  was  derived  and  verified  with  a  single 
element  compression  case.  With  the  beam  test  model,  the 
stiffness  reduction  technique  did  not  dramatically  affect 
the  material  with  fibers  oriented  ninety  degrees.  However, 
the  zero  degree  fiber  orientation  case  resulted  in  stress 
contour  changes  which  suggested  a  spreading  of  the  load 
through  an  element.  This  spreading  only  occurred  with  a 
stiffness  reduction  greater  than  ninety-five  percent.  At 
high  loads  and  stiffness  reductions,  the  model  breaks  down 
as  the  numerical  solutions  result  in  element  Instabilities. 
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